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1. Abstract 


Methods for iteratively determining the infrared optical constants for dust suspended 
in the Mars atmosphere are described. High quality spectra for wavenumbers from 
200 to 2000 cm -1 were obtained over a wide range of view angles by the Mariner 9 
spacecraft, when it observed a global Martian dust storm in 1971-2. In this research, 
theoretical spectra of the emergent intensity from Martian dust clouds are gener- 
ated using a 2-stream source-function radiative transfer code. The code computes 
the radiation field in a plane-parallel, vertically homogeneous, multiply scattering 
atmosphere. Calculated intensity spectra are compared with the actual spacecraft 
data to iteratively retrieve the optical properties and opacity of the dust, as well 
as the surface temperature of Mars at the time and location of each measurement. 
Many different particle size distributions are investigated to determine the best fit 
to the data. The particles are assumed spherical and the temperature profile was 
obtained from the CO 2 band shape. Given a reasonable initial guess for the in- 
dices of refraction, the searches converge in a well-behaved fashion, producing a fit 
with error of less than 1.2 K (rms) to the observed brightness spectra. The particle 
size distribution corresponding to the best fit was a lognormal distribution with 
a mean particle radius, r m = 0.66 /im, and variance, a 2 = 0.412 (r e //= 1.85 /im, 
v e ff=. 51), in close agreement with the size distribution found to be the best fit in 
the visible wavelengths in recent studies. The optical properties and the associated 
single scattering properties are shown to be a significant improvement over those 
used in existing models by demonstrating the effects of the new properties both on 
heating rates of the Mars atmosphere and in example spectral retrieval of surface 
characteristics from emission spectra. 
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2. Methods 


The goal of this work was to invert the Mariner 9 spacecraft data directly to deter- 
mine the optical properties of the dust suspended in the Mars atmosphere. 

To try to match any given Mariner 9 observation, we must know or estimate 
several things about the environmental conditions at the time of each measurement. 
The Mie and radiative transfer routines require these to be known in order to calcu- 
late brightness temperatures from given indices of refraction. Specifically, we must 
know the Martian surface temperature, the atmospheric temperature and pressure 
profiles, the dust optical depth and vertical mixing, and the dust particle shape 
and size distribution. None of these was independently measured at the time that 
the IRIS data were taken, so all must be inferred or estimated from the data, the 
literature or later spacecraft observations. 

2.1 Overview of Solution Technique 

Figure 2.3 illustrates the procedure used to iteratively determine the indices of 
refraction of the dust. Briefly, an initial guess for the imaginary index of refraction 
as a function of wavelength is made. A Kramers-Kronig routine is used to calculate 
the corresponding real indices of refraction. A lognormal particle size distribution 
is assumed for each data set and single scattering albedos are calculated using the 
Mie code. A multilayer radiative transfer code is called to calculate the theoreti- 
cal brightness temperature of the emergent radiation for each orbit. Atmospheric 
temperature profiles are taken from previously inverted C02 bands in the Mariner 
9 spectra, and we assume that the mixing ratio of the dust is uniform. Particle 
mean radius and the variance of the distribution are adjusted by finding the pa- 
rameters corresponding to the closest fit in the 550 cm -1 band of the data for all 
orbits. The temperature of the Martian surface and the optical depth for each orbit 
are adjusted by fitting the maximum transmission point (the wavenumber for which 
there is minimum absorption nor refraction) and the point of maximum absorption 
respectively. After the spectrum-specific parameters have been adjusted by fitting 
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appropriate regions of the spectra generated by the initial guess of optical properties 
with the spacecraft data, the optical properties themselves are then adjusted until 
the brightness temperature matches the data at all wavelengths and for all orbits. 

The new optical properties axe then made self-consistent using the Kramers- 
Kronig relations, and the entire process of adjustment is repeated until a set of 
self-consistent optical properties, particle size distribution, surface temperatures, 
temperature profiles, and optical depths is obtained. 

In this chapter we describe in detail the theory and computational procedures 
that allow us to do what was just described. 

2.2 Physical Model 

This section explains the relationship between the optical properties of the 
Mars dust, the physical parameters associated with the environment, and the inten- 
sity of thermal emission observed by an orbiting instrument. 

We begin with an introduction to the basic problem. As a radiation field 
propagates through an atmosphere, it is modified by various physical processes. In 
this work, we seek to isolate and study the interaction of atmospheric dust with 
the thermal emission in order to infer certain physical properties of the surface and 
dust particles. The molecular composition and the presence of aerosols such as dust 
or ice define the absorption and scattering properties of the atmosphere, which, 
combined with its thermal structure, determine the rate at which an atmosphere 
emits radiation at each wavelength in the thermal infrared. 

Although our method of inverting the data is unique and new, we use previ- 
ously well-established physical and computational methods to determine the radiar 
tive intensity emitted from the top of a plane parallel multiply scattering atmosphere 
[Toon et al. , 1977, Toon et al. , 1989, Goody and Yung, 1994], for the Mie single 
scattering theory which relates the optical properties of a spherical particle to its 
scattering properties, and for the relationship between the real and imaginary parts 
of the refractive index of a material [Hansen and Travis, 1974, Van de Hulst, 1971, 
Bohren and Huffman, 1983]. Here we review the physics of atmospheric radiative 
transfer and models used, and then outline our method of calling on these models 
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iteratively. 

Several process modify the radiation field as it passes through a dust laden 
atmosphere along a particular line of sight of an observer in orbit. First is the 
emission from the surface of the planet, which we assume to be radiating as a 
blackbody. This radiation is attenuated as it passes through the atmosphere due 
to extinction (absorption plus scattering out of the path) by particles in the path. 
Second is the emission of the dust particles themselves which are in the line of 
sight between the spacecraft and the ground. These particles are radiating with an 
emissivity that varies with wavelength and this radiation is also attenuated by the 
fraction of the atmosphere between the particle and the spacecraft. Lastly, for a 
scattering atmosphere, radiation from elsewhere in the atmosphere (other particles, 
other places on the ground, etc.) is scattered into the line of sight by particles in 
the path. The spacecraft measures the sum of all these intensities to obtain a total 
intensity as a function of wavenumber at that particular view angle. 

Here we shall briefly examine the general solution of the governing equation 
term by term to illustrate the above mentioned processes. 

Consider a plane parallel atmosphere. Here we define for convenience the 
vertical optical depth, r, which can be thought of as the optical cross section of a 
unit area of height d z, 

dr = —N 0 r e /fdz, (2.1) 


where N a is the total number of dust particles in the column and r e / / is the effective 
optical crossection of the particles. Note that the sign in this definition means that 
the optical depth goes from zero at the top of the atmosphere, to r at the bottom, 
with z defined as shown. The equation of transfer for an arbitrary, monochromatic 
field of radiation is then given by 


^ d t v 


-(1 - zu)B u (T) 

L I- 1 


( 2 . 2 ) 
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where © is the angle between the directions of incidence, (p, 0), and scattering 
(p',9 1 ), p is the phase function for single scattering, defined as the angular distribu- 
tion of radiation scattered once through the angle 0, and zu is the albedo for single 
scattering, defined as the ratio of radiant power scattered in all directions to that 
extinguished, or 

r4ir do/ 

W = Jo ^ ( 2 - 3 ) 
Physically, vo can be interpreted as the probability that a photon intercepted by a 
particle is scattered by that particle. 

The subscript v will heretofore be dropped for ease of notation, but it should 
remain understood that all parameters, in general, are functions of u, We assume 
azimuthal symmetry, implying that there are no horizontal gradients in atmospheric 
properties, so the dependence on <j> is eliminated. 

Our formal solution for the upward intensity emitted in the direction (p) from 
the top of the atmosphere is, term by term: 


= 

(1) /(r„/z)e-^ 

(2) + f (1 — t37)e _T ^ /1 S(r , )d— 

Jo 

( 3 ) + 27 T f f T - 

Jo J— 1 /i 

( 4 ) + ( 2 . 4 ) 

47 r Jo /i 

where the numbered terms represent: 

1. radiation originating from the planet’s surface, attenuated by extinction as it 
passes through the dust column of vertical optical depth, r* in the direction 
of p. 

2. radiation originating from a volume of emitting dust particles at a level r' , 
attenuated by the column of dust above the volume 


3. radiation from all directions //, scattered one or more times into the direction 
p by a volume of particles at a level r', attenuated by the column of dust 



above the volume 


4. radiation from the sun that has penetrated to depth r 1 and has been scattered 
back into the line of sight 


Note that the Planck function B is written B(r) to emphasize that tempera- 
ture is a function of optical depth alone (not x or y) in plane-parallel atmospheres. 
In this work, we neglect the last term (the solar term) because the solar flux in 
infrared wavelengths is negligible compared to the planetary radiation. 

Using approximations for the scattered intensity field and the scattering phase 
function in term 3, and employing the appropriate boundary conditions, the above 
equation is numerically integrated to calculate the upward intensity leaving the 
top of the atmosphere, given a particular set of surface and atmospheric tempera- 
tures, and given information about the absorbing and scattering properties of the 
atmosphere. Further, a vertically non-homogeneous atmosphere can be modelled 
by stacking individual layers and using the upward intensity from top of a lower 
layer as the boundary condition for the intensity entering the bottom of the layer 
above it. This stacking of layers will be discussed later in the section describing the 
multilayer code. 

The particles of dust in the atmosphere are of a diameter comparable to the 
wavelenths of radiation with which they interact. These particles are large enough 
to contain an optical ray, and we must treat all effects of diffraction, refraction, 
reflection, and phase shifts of the radiation as it passes through and around each 
particle. To relate the dust’s intrinsic material properties to its radiative scattering 
properties, we use the method of Mie [1908], which is an exact solution of the scalar 
wave equation in spherical coordinates 


1 d^(rip) 
r dr 2 


1 d . . £%>. 1 

t 2 sin 6 d6 Sm dd r 2 sin 2 6 d<ffi 


+ k 2/ ip = 0, 


(2.5) 


where the propagation constant, k is given by 
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Linear combinations of the resulting vector fields satisfy Maxwell’s equations. For 
a detailed treatment of the Mie theory, we refer the reader to any radiative transfer 
text [eg, Van de Hulst, 1957]. The Mie theory requires only that the particles be 
spherical, which is a good, but not perfect, assumption for the dust aerosols. 

The scattering properties we wish to calculate using the Mie theory are defined 
by the efficiency of the particles to scatter and absorb radiation. A particle of radius 
r has a geometric crossection, a g of 7rr 2 , and a wavelength-dependent ’’crossections” 
for extinction, absorption, and scattering, a ext , (Tabs, and a sct , respectively. These 
optical crossections can be thought of physically as the crossection of a particle that 
an incident wave of wavelength A sees, which may or may not be the same as the 
geometric crossection. We then define the efficiencies for extinction, absorption, and 
scattering as 


Qext = 
Qabs = 

Qsct = 


(Text 

°9 

(Tabs 


(Tsct 


(2.7) 


noting that extinction is the sum of absorption and scattering, so that 


(Text — ®a bs "t & Set 
Qext = Qabs "h Qsct 

According to Kirchoff’s law, the absorption crossection and efficiency factor of a 
particle are equal to its emission crossection and efficiency. 

The value of Q ext is typically on the order of unity, but can range from very 
small (-C 1) if the radius of the particle is much smaller than the wavelength of 
incident radiation (Rayleigh regime), to up to 4 or more for non-absorbing particles 
when of large radius (geometric limit). Since the scattering and absorbing properties 
of a given material are so sensitive to the ratio of particle size to incident wavelength, 
a study of the behavior of Q as a function of this ratio sheds valuable insight into 
the choice of particle size distribution. 
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It is useful to define a size parameter, 


for use in discussing the behavior of Q e xt- The particles in a dusty planetary atmo- 
sphere are typically not of one uniform radius. Instead, there is a distribution of 
particle sizes with a probability function n(r), with particles ranging in size in the 
case of Mars over several orders of magnitude (.001 /xm j r j 100 fim). 

The single scattering albedo, defined earlier as the ratio of radiant energy 
scattered to that extinguished, can be expressed in terms of the scattering and 
extinction coefficients: 


Finally, in discussing the scattering properties of a material, we would like to 
know something about the mean angle of scattering of a particle. This is represented 
by the asymmetry parameter, 


g = f (cos 0)p(cos<9)^, 
J 4ir 47T 


where p(cos 6) is the phase function we mentioned earlier. The phase function is 


normalized and given by 



( 2 . 11 ) 


and describes the probability of a particle to scatter radiation into the direction 
cos#. The asymmetry parameter, then, can have values between -1 and 1, where 
a value of +1 indicates that all radiation is scattered into the forward direction, 
a value of 0 indicates that radiation is scattered isotropically, and a value of -1 
represents complete backscattering. Typical values for dust particles such as those 
found in the Mars atmosphere are 0.1 < g < 0.8, with longer wavelengths being 
scattered more isotherm ally and shorter wavelenths being more forward scattered. 

The relationship between scattering efficiency and total optical depth is simple 
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and conducive to a physical understanding of these parameters. We may write 


T — N 0 GgQext 


( 2 . 12 ) 


or we can consider the mass of a unit area column of atmosphere 


roo 4 

M = I (p-Trr 3 )Nn(r)dr 
Jo 3 

= -pNreffCTg 


(2.13) 


where N is the total number particles of all sizes and p is the volumetric mass of a 
particle. We can then relate the optical depth and mass 


^ 3MQ ex t 

4pr eff 


(2.14) 


It is evident from equation 2.12 that since the number of particles, N, and the 
geometric crossection, <r g , are independent of wavenumber, we can find the optical 
depth at any wavenumber if we know the extinction efficiency at that wavenumber 
by a simple scaling using Q ext and r at a reference wavelength, or 


t(v) 


t(u = V ref ) 

Qex&iy = ^ref ) 


Qext (y) 


(2.15) 


We make frequent use of this convenient scaling to compute r quickly at all wave- 
lengths once it is known at a reference wavelength. 

The temperatures of both the atmosphere and the surface greatly affect the 
intensity of the radiation being emitted from the top of the atmosphere. The rela- 
tionship between these temperatures and the upward intensity is straightforward and 
intuitive, although the mathematical relationship is not linear or simple. Physically 
speaking, a warmer surface temperature will cause the brightness temperatures at 
all wavenumbers to increase, but there will be a more pronounced effect at wavenum- 
bers for which the indices of refraction of the dust are small. At wavenumbers for 
which the dust is optically thick due to high indices of refraction, the surface temper- 
ature has less effect because the surface is obscured by the intervening atmosphere. 
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Conversely, the temperature profile in the atmosphere has the greatest effect inside 
the dust absorption bands in the spectrum, with the temperature at the top of the 
dust layer effecting the temperatures at the cores of the bands and temperatures 
farther down affecting the wings of the bands. Therefore, the thermal structure of 
the atmosphere and the temperature at which the surface is radiating are important 
parameters in any radiative transfer model. The temperatures manifest themselves 
in Equation 2.4 in the Planck functions that are used to describe the intensity of 
radiation coming from the surface, as well as in the Planck functions for the emitting 
particles suspended in the atmosphere. 

One additional question we would like to answer concerns the overall heating 
or cooling effect of the suspended dust. A simple way to look at this to first order is 
to simply integrate the upward and downward intensities, I + and I~, which we have 
computed for each wavelength at specific emission angles, so that we may obtain a 
hemispherical mean flux quantity as a function of altitude, or 

F + (z ) = / 7T I + du 

Jo 

roc 

F~(z) = / 7rJ~di/. 

Jo 

Qualitatively, we would expect the upward flux to decrease with height from a 
value of <tT 4 at the surface to a somewhat smaller value at the top of the atmosphere. 
The downward fluxes also decrease but are much smaller in magnitude since they 
are due solely to downward scattering of radiation. These downward fluxes are zero 
at the top of the atmosphere and some small value at the surface. 

In summary, the set of parameters, Qext,m,g, constitute the scattering prop- 
erties of a material. They are all dependent upon wavelength, and also vary strongly 
with particle size, shape, and composition. The Mie theory provides an exact re- 
lationship between the optical properties of a particle and the particle’s absorption 
and extinction efficiencies. For a distribution of particles, an average set of efficien- 
cies can be found by integrating over all particle radii. The scattering properties 
therefore appear in the radiative transfer equation in the form of ru, e (the emissivity, 
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which, by Kirchoff’s Law, is given by 


e = (1 — w), (2.16) 

and p(cosO). The optical properties, rii and n r , plus information about particle 
size, shape, and distribution are contained within these scattering properties. The 
wavelength-dependent optical depth of the atmosphere appears directly in the equa- 
tion and can be computed at every wavelength using the scattering efficiency and one 
reference optical depth. Finally, the temperatures of the surface and atmosphere en- 
ter the radiative transfer equations inside the Planck functions for the surface (used 
in term (1) of Equation 2.4) and the emitting dust particles (term (2)). Overall 
heating effects of the dust can be examined by integration of upward and downward 
intensities over all wavelengths. 

The reader should now have a basic idea of the main parameters effecting our 
analysis of the interaction between suspended Mars dust and the thermal infrared 
radiation field of the planet. This brings us to the detailed discussion of the model 
and the methods for determining the parameters and properties of interest. 

2.3 Computational Models 

We use an iterative procedure to determine the indices of refraction of the 
suspended Mars dust, optical depths of the atmosphere for each orbit, surface tem- 
peratures for each orbit, and the particle size distribution of the dust. A detailed 
flow diagram showing the various components of the procedure is shown in Fig- 
ure 2.4. Atmospheric temperatures, surface emissivity, particle shape, and vertical 
mixing ratio are assumed to be known, and these assumptions will be discussed at 
length in following detailed descriptions of the model and procedures. 

The major subroutines of the model are the Mie single-scattering code, which 
calls on the Mie theory for each particle size radius then computes average scatter- 
ing properties for the dusty atmosphere; the Kramers-Kronig code used to compute 
physically consistent real indices of refraction for a set of imaginary indices; and the 
multilayer radiative transfer code which calculates the upward intensity of radiation 
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after it has passed through a vertically inhomogeneous, multiply scattering atmo- 
sphere by stacking individual layers of atmosphere and solving the radiative transfer 
equation for all layers simultaneously as a matrix of coupled equations. 

Procedures for using the above models consist of two broad categories of fits. 
First, the particle size distribution, surface temperature, and optical depth are sys- 
tematically varied, holding the indices of refraction constant at their initial guess, 
until the theoretical brightness temperature spectra generated by these parameters 
fit the data at specific wavenumbers. We refer to this process as the “adjustment” 
of particle size distribution, optical depth, and surface temperature. 

The second broad category of fits involves holding all the parameters which 
have just been adjusted fixed and searching at each wavelength for indices of refrac- 
tion that correspond to the best match of brightness temperature to the spacecraft 
data. This general process of fitting optical processes is sometimes referred to here as 
the “search”. Once new sets of optical properties are found, the real and imaginary 
indices must be re-correlated to each other using Kramers Kronig, which comprises 
the next initial guess of optical properties to be used by the adjustment routine. 

The above procedure is iterated until the optical properties which come out 
of the final search routine match those which were input to the adjustment routine. 

Now we will describe each computational model and procedure. 

2.3.1 Mie Code 

We have developed a driver program, called xmie, which calls a standard Mie 
code written at NASA Ames Research Center for every particle radius in a given 
particle size distribution at a given wavenumber. Subroutine xmie then integrates 
the scattering properties over all radii and divides by the appropriate normalization 
constant to obtain the average scattering properties for the layer of dust in question. 
It is capable of employing two different types of particle size distributions: the 
lognormal particle size distribution, which is the one we use in this thesis, and 
the modified gamma distribution [Diermendjian, 1969], which we can compute for 
comparison to previous work [Toon et al. , 1977, Clancy et al. , 1995, Ockert-Bell 
et al. , 1997, Pollack et al. , 1991]. 
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2.3. 1.1 Particle Size Distributions 

The logarithmic normal particle size distribution is a useful model for many 
polydisperse populations in nature [Campbell, 1995]. It is the model chosen to 
represent the distribution of dust particles in the atmosphere of Mars in this work. 

We assume the dust to be comprised of lognormally distributed particles for 
which the probability density function is 


n{r)dr 


1 

/o M/2 

ro^Z'Kyi 1 


' (In r-n) r 

2a 2 


dlnr 


(2.17) 


where // and o are the mean and variance of ln(r) respectively. It is normalized such 
that 

/•OO 

N(r) = / n(r)dr = 1.0 (2.18) 

Jo 

The distribution can be uniquely defined by the two parameters, fi and <r 2 , which 
describe logarithmic quantities and are therefore dimensionless. The notation of a 
for the variance of particle size distribution is unfortunately similar to the notation 
we’ve chosen for the geometric and optical crossections of the particles, but they 
should not be confused. 

Given the average geometric particle size, d = / 0 °° 7rr 2 n(r)dr, the average 
extinction and scattering properties of the atmosphere can then be computed: 


bar Q ^ 


XT «r 2 Q»ct(r)n(r)dr 

XT *rM r ) d r 


| barQ^ = 


XT’ 

XT^Mrldr 


barQ^yx — 


I , 00 *r 2 Q„ t (r)n(r)dr 

IT ,rr2n ( r ) dr 


(2.19) 


barg = 


XT 7rr 2 ff (r)n(r)dr 
X) 00 7rr 2 n(r)dr 


Similarly, 


( 2 . 20 ) 
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and 


w = 


Qsct 

Qext 


( 2 . 21 ) 


Figure 2.1 illustrates several lognormal distributions with different values of // 
and a. 

For verification of our integration of particle size distributions, we turning our 
attention to the moments of the probability function, 


r 


roc 

= / r n n(r)dr, 
Jo 


( 2 . 22 ) 


an analytical integration gives 


f n = [r] n e (r3n2/2 


(2.23) 


where n = 0, 1, 2..., we have the following: 

n = 0 f° = 1 
n = l f 1 = [r]e* 2 / 2 
n = 2 f 2 = [r]^ 
n = 3 f 3 = Me 9 *’/ 2 

These yield median and mean values 

| [r] = e^ median radius 
[s] = 7T [r] 2 median area 
[v] = 27r[r] 3 median volume 


(2.24) 


(2.25) 


(r) = f mean radius 

< (s) = 7 rf 2 mean radius (2.26) 

(v) = 27 rf 3 mean volume 

The code used in this work computes the particle size distribution, given jj, 
and a in subroutine xmie. At each particle radius, the effective radius and vari- 
ance are computed, along with the single scattering properties returned by xmie, by 
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Figure 2.1: Four particle size distributions are shown here: The dotted 
line is r 0 = .33, a = .8367, with r e ff and v e // matching Clancy et al. [1995]; 
The dot-dashed line is r 0 = 1.23, a = .5658, with r e ff and v e // matching 
Toon et al. [1977]; The solid line is r 0 = 1 ,<t = .1 
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integrating over 0 < r < oo. A gaussian quadrature routine, as described in Numer- 
ical Recipes [Press et al. 1996], pp. 121-125, is used to calculate the abcissas and 
weights at equal spacing in logspace for the integration. The analytical moments 
can be computed and compared with the integration to ensure accuracy. Using 
100 gauss points, the computed integration agrees with the analytical moments to 
within 10 -7 . 

It is common to see in the literature references to the particle “effective radius” , 
which is defined as the cross sectional area weighted radius, and a corresponding 
“effective variance” 

' _ _ XT rirr 7 n(r)dr _ f 3 

e f f 7T r 2 n(r)dr 

< . (2.27) 

_ / n °° r*r a n(r)dr 
V eff ~ / 0 °° wr a n(r)dr 

For comparison to previous research, the model can employ a modified gamma 
distribution, given by 

n(r) = Cr a <T*^ v (2.28) 

where r m is the mode radius of the distribution, a and 7 are positive constants, and 
C is the normalization constant to ensure / 0 °° n(r)dr = 1. 

The effective radius and effective variance, r^g and v e g, can be related analyt- 
ically to the lognormal mean radius and variance [Forget, 1996], r 0 and cr, 


r 0 = r e ff/(l + Weff/ 5 ^ 

a 2 = ln(l+Ueff), 


and to the modified gamma parameters, r m , a, and 7, by 


/'TM/ 7 r ( S ? 1 ) 
reff rm a T(^) 


Veff = 


r(«±5)r(a±3) 


r( s fr 


P^-i 


(2.29) 


(2.30) 


where T is the Gamma function T(x) = / 0 °° t x x e t dt. 

Although the lognormal and modified gamma distributions might look very 
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different from each other for the same values of r e g and Wes, Hansen and Travis 
[1974] showed that the average scattering properties depend very little on these 
differences and depend in most cases only on and v&. 

Ideally we would like to find a particle size distribution that matches observed 
data for all wavelengths in the electromagnetic spectrum. Past studies have had dif- 
ficulty finding particle size distributions that would match data both in the infrared 
and in the solar wavelenths. We seek in our model to find solutions for the optical 
properties using particle size distributions found by previous researchers to fit data 
in various wavelength regions [Toon et al. , 1977, Clancy et al. , 1991, Drossart et 
al. , 1991, Pollack et al. , 1995, Ockert-Bell et al. , 1997], and also to iteratively 
converge on a best fit of our own for comparison to these previously obtained parti- 
cle size distributions. We show results for cases in which particle size and variance 
were held fixed and all other free parameters adjusted and fit to the data, as well as 
results for the case when we let our model find its own best particle size distribution. 

2.3.2 Multilayer Radiative Transfer Code 

We use the word “multilayer” to indicate that our radiative transfer code 
is capable of approximating a vertically inhomogeneous atmosphere by stacking 
a number of single, homogeneous layers atop one another and solving the set of 
coupled equations to calculate the emergent intensity at the top of the stack. We 
will briefly sketch out the problem here. The atmosphere is divided into n layers 
(n = 4 here), numbered from top to bottom, j = 1,4. Atmospheric temperatures 
and optical depths must be specified at each boundary, so they are numbered from 
1 to n + 1 . Each layer is comprised of particles of particle size distribution rij (r) . We 
can write for each layer four equations, one each for J~, /+, I~, /+, and each layer 
has 6 unknowns. There are two equations of continuity between each layer, and 
two boundary conditions (one at the surface and one at the top of the top layer), 
therefore we have a total of 4*n+4*2 equations and the same number of unknowns, 
and we can solve directly for the intensity anywhere in the field. 

The model is capable of considering vertical inhomogeneities in particle size 
distribution, atmospheric temperature gradient, and optical depth, but is not ca- 
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pable of modelling any inhomogeneities in the basic dust particle composition (ie 
indices of refraction). A detailed study of the effects of varying particle distribution 
spatially and temporally would be of great interest, but is beyond the scope of this 
thesis. We assume for the purposes of our study that the particle size distribution is 
spatially and temporally constant. That is, we assume that the relative abundance 
of large to small particles is the same at higher altitudes as it is at lower altitudes, 
and that this distribution does not change throughout the duration of the duststorm. 
Both of these assumptions are contrary to physical intuition. Spatially, one might 
imagine that it’s very probable for larger particles to dominate the distribution at 
lower altitudes, and for smaller particles to dominate at higher altitudes, due to the 
tendency for particles to assume a vertical density profile according to the ratio of 
their potential energy to the thermal energy KT (where K is Boltzmann’s constant 
and T is the local equilibrium temperature). For this natural vertical particle size 
profile not to occur, there must be enough turbulent mixing as to keep the parti- 
cles uniformly distributed. We employ the the conclusion by Toon et al. [1977] and 
Conrath et al. [1975] that Stokes Cunningham fallout is not the dominant mecha- 
nism leading to the dust storm clearing and that strong vertical mixing persisted 
throughout the duration of the storm. This assumption, although supported by 
some studies, is controversial [Haberle, private communication], and merits further 
study. 

Although the size distribution of particles is assumed to be constant with 
height, the density of the particles is not. We assume a constant ratio of dust mass 
to atmospheric mass, so that the particles become less dense with height. 

Table 2.1 gives list of physical parameters that must be specified to generate 
a synthetic spectrum. 

Several assumptions and approximations are necessary to make the solution 
of the multilayer radiative transfer equation tractable. Briefly, the most important 
are 

1. the representation of the radiation intensity field inside term (3) of Equa- 
tion 2.4 as two streams of radiation, one in the positive (upward) direction 
and one in the negative (downward) direction, 
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Table 2.1: Parameters used in solution of multilayer radiative transfer 
equation 


Symbol 

Parameter 

Specified by 

nlayer 

nspec 

ndat 

rii(ndat) 

n r (ndat) 

T 8UT t (nspec) 
Ti(nlayer + 1 , nspec) 
tau^(nlayer + 1, nspec) 
r m (nlayer, nspec) 
a (nlayer, nspec) 
e 8UT f (nspec) 

Number of layers 
Number of spectra 
Number of datapoints 
Imaginary indices of refraction 
Real indices of refraction 
Surface Temperatures 
Atmospheric Temperatures 
Ref. Optical Depth, A = 9pm 
Particle mean radius 
Particle size dispersion 
Surface emissivity 

Choice 

Choice 

Smoothed Mariner 9 data 
Iterative bisection 
Kramers- Kronig 
Iterative bisection 
C0 2 band inversion 
Iterative bisection 
Step search 
Step search 
Assumed black=1.0 


2. use of just the first terms of the legendre polynomial expansion of the phase 
function, p(cos0), also inside term (3), 

3. use of just the first term of the Taylor expansion of B(r) in determining the 
Planck function in terms (1) and (2) of Equation 2.4, namely 

B(t) = B 0 + B 1 t (2.31) 

4. use of the approximate atmospheric temperature profiles inverted from the 
C0 2 band using the assumption that dust was not absorbing in that region 

5. calculations and searches performed only at wavenumbers outside the major 
C0 2 absorption bands 

6. the assumption that incoming solar radiation is negligible in the infrared, 
compared radiation from the Martian surface 

7. the assumptions that the dust is well mixed from the ground up to a height 
where the atmospheric pressure is 0.1 mbar and that the mixing ratio of dust 
to atmosphere remains constant to that height, where the dust ends abruptly 

8. the approximation that dust particles behave as spherical particles. 
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9. the approximation of surface emissivity of unity rather than having a spectrum 

of its own [Christensen et al. , 1997] 

The multilayer code returns the upward radiation in the form of brightness 
temperature for comparison to the spacecraft data. 

2.3.3 Searching routines 

A number of searches are performed on various parameters to obtain the best 
fit to the data. Parameters for which searches are performed include the reference 
optical depths, the surface temperatures, the particle size distribution, and the 
optical properties themselves. In all cases except for the particle size distribution, 
the criterion for fitting is the square of the error between the computed and the 
actual brightness temperatures. For properties which are the same for all orbits 
(optical properties and particle size distribution), the square of the error is summed 
over all orbits. In all cases except for the optical properties, the error function, 
which we call x, is well behaved, and a simple bracketed bisection searching method 
is used, holding all other parameters besides the one sought constant and finding the 
value of the parameter of interest which generates the minimum error in brightness 
temperature. 

2. 3. 3.1 Particle Size Distribution 

The particle size distribution is the only search which is not completely auto- 
mated. Two separate loops stepping through a range of possible mean particle radii 
and variances are wrapped around the outside of the code, and the best particle 
size parameters are chosen after the rest of the searches have been carried out for 
each possible combination. This is a long and laborious process, which we often 
disable and instead pick a reasonable particle size distribution a priori based on the 
literature, and perform the subsequent searches on just that one set of particle size 
parameters. 

2. 3. 3. 2 Surface Temperature and Tau 

Surface temperature and reference optical depth are strongly coupled, so the 
process of searching on these parameters is a smaller iterative loop within the 
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larger “adjust” iteration. First, using the initial guess for optical properties at the 
wavenumber for which the brightness temperature is a maximum (discussed ear- 
lier as the wavenumber most sensitive to surface temperature), and using an initial 
guess for reference optical depth, holding all other parameters constant, the surface 
temperatures are found for each orbit which produce brightness temperatures that 
match the spacecraft data at that one maximum point in the curve. Then, using 
those new values for surface temperature and holding everything else constant, the 
reference optical depths are adjusted until the curves fit at the waveneumber in the 
core of the 9 //m dust band. This fit causes the curves at the maximum brightness 
temperature point not to match anymore, so another search on surface temperature 
is done, and so on until both points fit exactly for all orbits. This iteration always 
converges as long as the intial guess for optical properties has at least a rough dust- 
like absorption feature, so that the brightness temperatures are more sensitive to 
the surface temperature than to optical depth at the maximum brightness temper- 
ature point, and less sensitive to surface temperature than to optical depth at the 
wavenumber corresponding to the core of the dust absorption feature in the data. 

2. 3.3.3 Indices of Refraction 

The complex index of refraction is given by 

n( A) = n r ( A) — irii( A) (2.32) 

where n r is the usual refractive index which determines phase speed, and n t is related 
to the absorption coefficient 

k a =A'KTiij\. (2.33) 

The real and imaginary parts of the index are not independent, but connected by in- 
tegral relations. Given one set of indices over a sufficiently large band of frequencies, 
the other can be obtained by simple integration. For a derivation of these relations, 
please see Bohren and Huffman [1983], pp. 19-22. The Kramers-Kronig relations 
(also known as the dispersion relations) relate the real and imaginary parts of the 
complex index of refraction of a material, based on the assumptions of linearity 
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and causality in the interaction of a frequency dependent electric field with matter. 
We compute the real index of refraction of the dust particles at each wavenumber 
using information about the imaginary indices of refraction at all wavenumbers in 
the spectrum. The dispersion relations can be written in many different forms, and 
many different methods for solving the integral equations have been used [Warren, 
1984], 

The following subtractive-Kramers-Kronig (SKK) relation is chosen for our 
calculations since it is the preferred method if 7i< m is known only in a short spectral 
region [Warren, 1984, Ahrenkiel,1971, Bachrach and Brown, 1970, Hale and Querry, 
1973]. 

n r (Ao) - n r (Ai) + * P J q _ A 2)( A 2 _ A 2 ) d ( ln A )- ( 2 34 ) 

where P is the Cauchy Principle Value of the integral. The above can be rewritten 
as 


2 

(A 0 ) = n r (Ai) + - / A 2 n im (A) 
7T J 0 


n< 


[(A* -A?) (A2- Ag)J 


cf(lnA). (2.35) 


In our model, we provide the Kramers-Kronig code with a limited range of imaginary 
indices of refraction (thermal infrared, 5-40 //m), and give it a known value of the 
real index at a wavelength in the visible region of the spectrum (Ao=0.5 //m), for 
which the properties have been well established [Ockert-Bell et al. , 1997]. 

We use a modified version of the code written by Sinnot at NASA Ames 
Research Center [Pollack, private communication]. The subroutine numerically in- 
tegrates the integral in pieces on both sides of each singularity and sums them. It 
interpolates linearly if the wavelength is within the range for which data exist and 
extrapolates if not, according to the following: 


72j(A) ~ nii^smallest ) A^ forA ^ A smallest 
ft»(A) ~ ^t(Aiarpe«t) A for A ^ A largest 


(2.36) 


The subroutine takes as input an array of wavenumbers (in decreasing order), a 
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set of imaginary indices of refraction at the given wavenumbers, a value for the 
corresponding real index of refraction at one reference wavelength, the size of the 
wavelength array, and an integration resolution factor. It returns the array of real 
indices of refraction corresponding to the imaginary set it was given. Figure 2.2 
shows the accuracy of the Kramers-Kronig subroutine. The bottom plot is the set 
of imaginary indices, and the top compares the actual real indices with the SKK 
computed ones (upper curve). The rms error between the calculated and actual in- 
dices is 6.6 x 10 -2 . Removing the carbon dioxide band from the integration reduces 
the accuracy in the calculation, particular in wavenumber regions immediately ad- 
jacent to the 15 /im band. The error caused by linear inerpolation of this band in 
the calculations is shown. 

Assuming the temperatures, optical depth, and particle size distribution to 
be known, there can exist at every wavelength more than one complex index of 
refraction which will generate a given brightness temperature. If there are two 
solutions, one represents the scattering solution and the other the absorbing solution. 
Although not every complex index is physically feasible, this non-uniqueness makes 
an inversion of the data problematic. The phase space is not well suited the type of 
search necessary to find one value for ni which produces minimum error. A technique 
for simultaneously matching multiple orbital spectra to eliminate the ambiguity was 
developed [Toon, personal communication]. One unique set of imaginary and real 
indices of refraction can be found which generates brightness spectra that match 
multiple data sets simultaneously. 

The same type of bracketed bisection search is used for the optical properties. 
The code can either search only on rii, which can sometimes cause difficulties in 
non-absorbing regions where the brightness temperature is more sensitive to the 
real index than to the imaginary, or it can search on both n* and n r in parallel 
then combine them gradually to iteratively obtain a physically consistent set of 
properties. 
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Figure 2.2: Real indices of refraction for Montmorillonite 219b [Toon 

et al. 1977] calculated with Kramers Kronig by integrating over all 
wavenumbers (250-2000 cm -1 ) (dashed) and calculated without the C0 2 
band (dotted) compared with the actual values (solid). 
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Table 2.2: Test Parameters, (Ockert-Bell et al. [1997] particle size distri- 
bution) 


Optical Depth 
(reference) 
i/=1075 cm -1 

View 

Angle 

cos© 

Bottom Atm. 
Temp. 

(K) 

Top Atm. 
Temp. 

(K) 

Surface 

Temp. 

(K) 

Mean particle 
radius 
r 0 (//m) 

Var(ln r 0 ) 

1.286 

0.417 

214.177 

243.746 

266.000 

0.660 

0.642 

1.417 

0.904 

199.792 

259.488 

260.000 

0.660 

0.642 

1.024 

0.969 

196.730 

258.875 

263.000 

0.660 

0.642 

0.761 

0.387 

191.079 

246.290 

278.000 

0.660 

0.642 

0.442 

0.996 

179.595 

249.974 

275.000 

0.660 

0.642 

0.424 

0.982 

178.460 

255.760 

278.000 

0.660 

0.642 

0.405 

0.527 

174.716 

245.766 

281.000 

0.660 

0.642 

0.349 

0.505 

168.652 

241.217 

278.000 

0.660 

0.642 

0.349 

0.980 

168.652 

241.217 

278.000 

0.660 

0.642 

0.235 

0.974 

169.684 

240.538 

281.000 

0.660 

0.642 


2.4 Test Case 

To demonstrate our techniques, we have constructed ten brightness tempera- 
ture spectra using montmorillonite 219b, for which the optical properties are known, 
using different optical depths and surface temperatures, as shown in Table 2.2. These 
values were chosen in order to verify our model with previous work, and to test con- 
vergence of the various search techniques. The ten synthetic spectra generated using 
these parameters and montmorillonite 219b are treated as they were to be matched, 
then the searches on surface temperature, optical depth, and optical properties are 
carried out and the results compared to the known values. The search on particle 
size was not verified using this test case. 







3. Results and Discussion 


Given a reasonable initial guess for the indices of refraction, the searches converge 
in a well-behaved fashion, producing a fit with error of less than 1.2 Kelvin (rms) 
to the observed brightness spectra. The particle size distribution corresponding to 
the best fit was a lognormal distribution with a mean particle radius, r m = 0.66, 
and variance, a 2 = 0.412 (r e //= 1.85, i/ e //=.51) (see Figure 3.7), in close agreement 
with the size distribution found to be the best fit in the visible wavelengths in recent 
studies by Pollack et al. , [1995] and Ockert-Bell et al. [1997], Figure 3.1 shows the 
final fit for all spectra (solid line) to the data (points) for each orbital dataset, 
beginning with the earliest (dustiest) orbit 008 in the upper left and ending with 
the latest (clear) orbit 206 in the lower right comer. The optical properties which 
were found to best match the data are shown in the lower plot of Figures 3.2- 3.4, 
with the dashed lines representing the initial input to the search, and the solid lines 
representing the converged result. The orbital parameters which correspond to these 
optical properties and which were also found iteratively during this search are shown 
in Table 3.1 and in Figure 3.5 by solid lines. Resulting scattering properties and 
opacities for one of the dustier orbits are shown in Figure 3.6. Since we have assumed 
a spatially and temporally constant particle size distribution and composition, these 
values of single scattering albedo and asymmetry factor represent all orbits, but the 
opacities scale for each orbit according to the density of the dust at the particular 
time of measurement. Again, the values found by Toon et al. . are shown for 
comparison. 

3.0.1 Uniqueness 

The problem examined here is as follows: given a measurement or series of 
measurements of thermal radiation emitted by an atmosphere, the intensity and 
spectral distribution of which depend on the state of the atmosphere in a known 
way, deduce the best estimate of the properties of the atmosphere. There are two 
distinct aspects to this problem which are not always clearly separated; they may 
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Figure 3.1: Fit produced by the new optical properties shown in Fig- 
ure 3.2 to brightness temperatures measured for orbits ranging from 
dusty (r(9/zm) = 1.5) (upper left) to clear (r = 0.18). The observations 
are represented by points and the fits are solid lines. 
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Figure 3.2: Upper plot shows initial (dashed line) and final fit (solid line) 
to the Mariner 9 spectra (points) for orbit 8. Lower plot gives the initial 
(dashed lines) and final (solid and dotted lines) optical properties which 
produced the temperatures in the upper plot. 
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Figure 3.3: Upper plot shows initial (dashed line) and final fit (solid line) 
to the Mariner 9 spectra (points) for orbit 80. Lower plot gives the initial 
(dashed lines) and final (solid and dotted lines) optical properties which 
produced the temperatures in the upper plot. 
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Figure 3.4: Upper plot shows initial (dashed line) and final fit (solid line) 
to the Mariner 9 spectra (points) for orbit 8. Lower plot gives the initial 
(dashed lines) and final (solid and dotted lines) optical properties which 
produced the temperatures in the upper plot. 
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Figure 3.5: The solid lines in the upper plot show the surface temperature 
(higher solid line) and the atmospheric temperature near the surface 
(lower solid line) during the decay of the dust storm. For comparison, 
the assumption made by Toon et al. is illustrated by the dashed line, and a 
very crude verification of trends obtained from the LMD database [Forget 
et al. 1998] is given by the dotted lines. The lower plot compares the new 
inverted opticity at a wavelength of 9 ^im with the opacity assumed by 
Toon et al. , [1977]. 
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Figure 3.6: Final computed optical properties (top), with new properties 
shown by the solid lines, montmorillonite shown by the long dashes, and 
palagonite shown by the short dashes. Extinction efficiency, Qext, (middle 
left) single scattering albedo, m, (middle right) asymmetry parameter, g, 
(bottom left), and spectral opacity, rx, are compared for various studies. 
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be described as the ’inverse’ problem and the ’estimation’ problem. The inverse 
problem is the matter of inverting a known equation which expresses radiation as a 
function of the atmospheric state, so as to express the atmospheric state in terms 
of radiation. This is usually an ’ill-posed’ problem; i.e., it has no mathematically 
unique solution [Rodgers, 1972]. We therefore have an estimation problem, that 
is, to find the appropriate criteria which determine the best solution from all the 
possible ones which are consistent with the observations. 

We have dealt with this situation of non-uniqueness in two ways. First, we 
have made approximations to the equations of radiative transfer by discretizing 
continuous functions of altitude, wavenumber, and view angle, or by expanding 
functions in Taylor series or Legengre polynomials. Secondly, we have imposed 
constraints on a large number of the unknowns based on strong physical arguments 
which will be discussed in Section 3.0.2. 

Our resulting problem is one which consists of many more equations than 
unknowns, but the equations are very strongly coupled and non-orthogonal. Specif- 
ically, there exists an equation of radiative transfer which relates the our unknowns 
(the indices of refraction, the surface Temperatures, the particle size distribution, 
and the opacity) to the measured radiance quantity at each wavelength and for each 
orbit. There also exists at each wavelength one Kramers Kronig relation which re- 
lates the value of the real index of refraction of the material to the integral of the 
imaginary indices over all wavelengths. We have selected ten datasets of 85 wave- 
lengths over which to solve these equations. Thus, we have 10*85+85=935 total 
equations. Our remaining unknowns after discretization and physical constraints 
are 9 surface temperatures (one for each orbit), 9 reference optical depths, 85 imagi- 
nary indices of refraction, and 85 real indices of refraction, one mean particle radius, 
and one size distribution variance, for a total of 9+9+85+1+1= 105 total unknowns. 

Once the problem has been so reduced, we may treat it as well-posed in the 
immediate vicinity of the solution. We use the common technique of a least squares 
fit in the error between the computed and observed brightness temperatures. The 
105-dimensional phase space over which we must search is riddled with local minima. 
We search on each parameter individually by traversing one dimension of the phase 
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space at a time and finding the minimum in each dimension. The dimensions of 
phase space for ni, nr, r, and T s , are well behaved in the vicinity of the solution, 
with either one or two sharp minima which are easy to find. Provided we do not 
stray too far from the actual solution with our initial guess, the result of the iterative 
search scheme finds the appropriate solution. However, the scheme is not generally 
robust. That is to say, it is prone to either divergence or convergence to any number 
of mathematically feasible results, if given a bad initial set of parameters. This would 
not be a problem if all solutions besides the right one were physically unrealistic and 
could be ruled out on the basis of physical arguments. However, we have found that 
many physically feasible solutions exist. 

Fortunately, we know enough about the physics of the situation to ensure 
that we start with a guess which is close enough to the solution to avoid the above 
mentioned problems of divergence or wrong results. In addition, our data are noisy 
and are averaged over many measurements, which tends to damp out minima in the 
actual phase space, making the search more robust than it otherwise would be. 

To investigate the robustness of the search, we constructed a test case by 
using a known spectrum of optical properties to calculate brightness temperatures 
for various values of optical depth and surface/atmospheric temperatures and then 
attempted to invert the problem to obtain the values we started with. This artificial 
test case had an unfriendly phase space with sharp multiple minima which would 
trap the search in incorrect solutions and which even diverged from the correct 
solution if perturbed slightly. There seemed to be a finite number of incorrect 
solutions for the limited number of initial guesses we tried. That is, most initial 
guesses converged to one of three or general solutions. 

We then added noise to the test case and discovered that this greatly reduced 
the undesirable behavior of the search, and the model was capable of converging to 
the right answer given any reasonable initial guess. In fact, with noise comparable 
to that found in the actual data, it was possible to start with a flat gray guess for 
the imaginary indices of refraction (a very bad guess, indeed) and obtain the correct 
solution. 

Thus we are confident that while the solution we have obtained is not unique, 
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it is the correct one, since we know we are starting with an initial guess which is 
within the tolerance of proximity to the actual solution. Figures 3. 2-3.4 show in 
dotted lines our initial gray guess, and the final converged answer for the optical 
properties, along with the fits to the brightness temperature spectra. This solution 
is the one obtained by a large number of initial guesses, increasing our degree of 
certainty in the result. 

3.0. 2 Accuracy of the model and plausibility of the results 

Now that we are confident in the results of the code used to invert the spec- 
tra, we must convince ourselves that the assumptions made in order to solve the 
equations on which the code is based and the results which the code generates are 
justifiable and physically realistic. 

3.0. 2. 1 Particle Size distribution 

Of all the free parameters in this problem, the particle size distribution has 
proven to be the most difficult to invert directly. The minimum in the particle size 
distribution phase space is very broad and covers a large range of radii. We had 
originally intended to consider the particle size distribution parameters as totally 
free, and to try to invert them directly with the rest of the free parameters. However, 
limitations on the computational time combined with the large number of possible 
solutions in this dimension of phase space led us to conclude that this aspect of the 
problem was better suited to individual case studies rather than a sweeping broad 
search of all possible particle size distributions. 

For our case studies, we focused attention on particle size distributions which 
have been proposed in the literature [Hanel et al . , 1972, Hunt et al . , 1972, Conrath 
et al. , 1973, Toon et al. , 1977, Santee et al. 1993, Clancy et al. , 1995, Pollack et 
al. 1995, Forget et al. 1997, Ockert-Bell et al. 1997]. Some of the distributions which 
produced the best fit in our model are shown in Figure 3.7. We used single-mode 
lognormal distributions with equivalent effective radius and variances (see Chapter 2 
to the previously proposed modified gamma distributions. Hansen and Travis [1974] 
demonstrated that the effective radius and effective variance are the only important 
parameters of the size distribution, and that the radiative transfer equation is not 
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sensitive to differences between types of distributions (eg. the differences between 
lognormal and modified gamma type distributions) in the infrared wavelengths. 

All size distributions which were investigated represent physically plausible 
particles, although those proposed by Clancy et al. [1995] are probably too small, 
and those proposed by Toon et al. [1977] are probably too big. Multi-modal size 
distributions (representing mixtures of very different materials) were not considered. 

The particle size distribution which was ultimately chosen as the best was 
r e ff = 1.85, v e ]f = .51, also found by Ockert-Bell et al. [1997] and Pollack et 
al. [1995] to provide the best fit to spectra at visible wavelengths. Previous at- 
tempts to fit brightness temperature spectra of earth analog materials using this 
size distribution had failed to produce a suitable result in the infrared [Clancy et 
al. 1995], but we find it is possible to fit all brightness temperature spectra using our 
inverted optical properties and this particle size distribution. Thus, we now have 
a proposed size distribution which can fit data from the infrared to the ultraviolet 
extremes of the spectrum. This continuity of fit over all regions of the spectrum is 
a strong argument in favor of this particular distribution. 

Numerical integration over particle size in the calculation of the Mie scatter- 
ing parameters introduces only a small error in the calculations, since an accurate 
gaussian quadrature is used to determine weights and abcissas for the integration, 
and a high number of intervals is used 

Two other important assumptions regarding the particle size distribution were 
made in this analysis: l)particles were uniformly mixed throughout the atmosphere 
and 2) particle size distribution remained constant throughout the duration of the 
storm. The first assumption is probably an accurate one, based on calculations of 
wind speeds and eddy velocities [Toon et al. 1977, Haberle et al. 1996], showing that 
winds are strong enough to keep the dust well mixed up to altitudes of 50 km. Our 
results suggest that second assumption, however, may not be so accurate. Figure 3.1 
clearly shows an overprediction of brightness temperature on orbit 8 over the region 
of the spectrum most sensitive to particle size (450 cm -1 ), and an underprediction 
later in the storm on orbit 154, indicating that the particles might actually be bigger 
than assumed near the beginning of the storm and smaller than assumed towards 
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Figure 3.7: Lognormal particle size distributions with effective radii 

and effective variance corresponding to distributions used by Toon et 
al. [1977], Pollack et al. [1995], Clancy et al. [1995], Forget et al. [1997], and 
author]. 
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the later clearing stages. Nevertheless, the fits are still the best ever found in this 
region of the spectrum (see Section 3.0.4), and will only improve once variation in 
particle size with time is taken into consideration. 

3.0. 2.2 Optical Properties 

The only a priori assumption made regarding the optical properties is that 
the composition of the dust remains constant. That is, the density of the dust 
may change throughout the storm, but the dust has the same intrinsic material 
properties at the beginning of the storm as at the end. The data Figure 3.8 support 
this assumption, since the absorption and transmission features for all spectra occur 
at identical wavelengths, and all spectra are self-similar to first order. 

It is difficult to say how well the Kramers Kronig equations represent the actual 
correlation between the real and imaginary indices of refraction, since most mea- 
surements of optical properties in laboratories employ various theories, including 
Kramers Kronig, to obtain self-consistent real and imaginary indices of refraction. 
For example, all ’’actual” material properties employ Kramers Kronig relations, 
Habke theory, or dispersion analysis in some form [Rousch, 1994]. The Kramers 
Kronig equations are exact, derived from the Maxwell equations, so they repre- 
sent the best known relationship between the real and imaginary refractive indices. 
Again, the numerical integration of the relations introduces error near the region of 
the interpolated CO 2 band and near the endpoints of the region of integration, but 
the use of the subtractive forms of the relations reduces the effects of the truncation 
to the immediate vicinity of the endpoints. This has been discussed in greater detail 
in earlier sections and illustrated in Figure 2.2. 

During the iteration procedures, there are moments when the real and imag- 
inary indices of refraction are not Kramers Kronig consistent, but as the searches 
converge, this is no longer the case, and the resulting optical properties conform to 
the physical constraints imposed by the dispersion relations. 

The optical properties which result from the iteration scheme have already 
been shown in Figure 3.2. Figure 3.6 demonstrates by comparison (top left and 
right)to properties of the analog materials, montmorillonite and palagonite, that the 
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Figure 3.8: This figure shows the similarity of the Mariner 9 datasets 
chosen for ttiis analysis. In general, all 20,000 Mariner 9 spectra exhibit 
roughly the same behavior as shown here. The three plots show the 
successive smoothing employed, and also demonstrate the consistency of 
composition and material properties of the dust over time. 
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resulting indices of refraction realistic and exhibit the general absorption features 
expected in the region centered at approximately 1075 cm -1 and 550 cm -1 typical of 
silicate materials such as clays and basaltic weathering products. The problematic 
550 cm -1 double absorption peaks of the montmorillonite 219b imaginary index do 
not appear in our results. In addition, new absorption features appear in the small 
wavenumbers (< 500 cm -1 ) which have never been included in traditional silicate 
models of optical properties in the infrared. These features could be indicative 
of carbonates such as gypsum or calcite, which have similar properties at these 
wavelenghts [Rousch, private communication, Long et al. 1993]. Gypsum and calcite 
have been identified by Viking and Pathfinder to be possible components in the Mars 
surface soils and rocks, so it would not be suprising to find some evidence of these 
minerals in atmospheric dust spectra. It is somewhat surprising that the absorption 
does not decrease near the limit of small wavenumbers. This could be due in part 
to the high level of noise in the data at these wavenumbers and in part to possible 
inaccuracies in particle size and shape. Near this Rayleigh scattering limit, the 
assumptions of spherical particles begin to become less accurate [Pollack, 1982]. 

The real index of refraction behaves as expected, with a phase-like response to 
the imaginary peak and asymptoting to the known value for n T toward the visible 
end of the spectrum of 1.52 [Ockert-Bell, 1997]. 

3.0. 2. 3 Surface and Atmospheric Temperatures 

The surface temperatures on which the model converges are within 8 K of 
the approximate surface temperatures given by the brightness temperature at the 
maximum point (1321 cm— 1). The dashed line in Figure 3.5 shows the brightness 
temperature at the wavenumber for which the measured radiance is maximum, cor- 
responding to the wavenumber of minimum opacity of the dust. Toon et al. [1977] 
assumed this value to be the temperature of the suface, and if we assume that the 
atmosphere is cooler than the surface, this dashed line represents the coldest possi- 
ble surface temperature. We have not assumed that the dust is entirely transparent 
at any wavenumber, so the surface temperatures in our results (upper solid line in 
upper plot) are warmer than those assumed by Toon et al. . The lower solid line is 
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the value of the bottom temperature at the bottom of the atmosphere, just above 
the surface. Large temperature gradients between the surface and in the first few 
meters of the atmosphere have been observed by both Viking and Pathfinder , and 
we see this reflected in our results, particularly in the later clearer conditions, where 
the atmosphere above the surface is 40 Kelvin cooler than the surface itself. 

An additional check of the physical reality of the temperatures at the time 
and place of each measurement was performed using the European Mars Climate 
Database at http://www.lmd.jussieu.fr/mars.html [Forget et al . , 1998]. This database 
combines data from Mariner 9, Phobos, Viking, and Pathfinder, and uses theoreti- 
cal models for interpretation between datasets. The data are sparse and the models 
are coarse, but we obtained some rough verification of the trends in temperatures 
based on the database values for a dusty atmosphere (r = 2) during the northern 
winter season. The dotted lines in the upper plot of Figure 3.5 show the values 
returned from the database, using a constant value 2 for the opacity until orbit 154 
and abruptly switching to the “low dust” scenario after orbit 154. Clearly, this will 
return incorrect temperatures, since the opacity decreases smoothly from 1.5 to .2 
over the course of the storm. Nevertheless, the relative values of the temperatures 
show the effects of time of day and location of measurement, and the slopes of the 
lines, particularly in the dustier orbits, closely match our results. 

As discussed in Chapter 2, we do not consider the vertical temperature profile 
to be a free parameter. We use the temperatures inverted from the 667 cm -1 CO 2 
band by Pollack et al. [private communication], which do not take into consideration 
scattering or absorption by atmospheric dust. Other researchers [Santee and Crisp, 
1993, Lindner, 1993] have attempted to treat the dust and CO 2 simultaneously, but 
did not invert the particular spectra which we selected for our studies, and a detailed 
analysis of the combined C 02 /dust absorption is beyond the scope of this thesis. 
We therefore accept the error associated with the inverted temperatures, which is 
on the order of 5-10 Kelvin [Santee and Crisp, 1993]. Sensitivity studies involving 
these assumptions and errors are presented in Section 3.0.3. Since the radiative 
transfer equations are only very loosely coupled in wavelength, we are able to focus 
our attention on only those wavelengths outside the CO 2 band, and we make no 
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claims regarding the optical properties in that band. We treat that region as one 
in which there is no data, which we can do because the radiative transfer equations 
are essentially monochromatic, and the Kramers Kronig equations are not sensitive 
to wavelengths far away from the one being calculated. 

3.0. 2.4 Opacities 

The opacity of the dust is considered to be a completely free parameter, with no 
a priori assumptions. The optical depth at wavenumber 1075 cm -1 was iteratively 
determined for each orbit, and then opacities at all other wavelengths were computed 
using the extinction efficiency scaling ratio as discussed in Chapter 2. 

The bottom plot of Figure 3.5 gives a normalized comparison of optical depth 
decay rate with that assumed by Toon et al. [1977]. The solid line is the result 
of our calculations, as compared with the dotted line which has been assumed in 
most studies of the Mariner 9 data. These assumptions were based on knowledge of 
opacities in other parts of the spectrum (eg. Viking opacity measurements in the 
visible during 2 dust storms) and on data from other instruments aboard Mariner 
9. The trend in our results is clearly similar, but indicates that the dust probably 
persisted longer near the beginning of the storm until after orbit 60 when it began 
to decay. 


3.0. 2. 5 Other assumptions 

We have included the effects of a spectrally varying surface emissivity, to reflect 
the conclusions of Christensen et al. [1998] that the surface signatures are included in 
the observations of radiance and must be separated from atmospheric effects. This 
represents the physical situation better than assuming the surface to be a black or 
gray body, such as Toon, Pollack, Forget, and Clancy have assumed, but is only an 
approximation based on the spectra measured during clear atmospheric conditions 
late in the Mariner 9 mission. 

The multilayer radiative tranfer model was verified using simple tests for which 
results were known. 
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3.0. 3 Precision of calculations 

Brightness temperature fits for all ten orbits are shown in Figure 3.2 - Fig- 
ure 3.4. We find that the RMS error in brightness temperature is at or below the 
noise-equivalent AT for all orbits. We have done a systematic set of studies to de- 
termine the sensitivity of our model to the various parameters listed in Table 2.1. 
Table 3.2 shows the rms error for our baseline result (Table 2.1), with a breakdown 
of error for each orbit. Each parameter was increased and decreased by 10%, except 
in the case of the Temperatures, which were increased and decreased by 5 K. All 
other parameters were held constant in all cases. We see that the rms error increases 
in all cases as we move away from our baseline answer, indicating that we have in- 
deed found the best fit of these parameters for this given particle size distribution 
(reff = 1.85, Vfff = 0.51, and giving an idea of the error bars for the results. 

3.0. 4 Implications of the results 

There are two general categories of studies for which improved optical proper- 
ties are necessary. First, current spacecraft measurements from orbit must be able 
to remove atmospheric effects in order to study the spectral properties of the surface. 
Current spectral models for the properties of the atmospheric dust in the infrared are 
inaccurate, as we will demonstrate. The atmosphere is always dusty, and therefore 
there will always be a need for accurate characterization of the spectral properties 
of the dust. Secondly, any measurements involving heating rate calculations of the 
surface or atmosphere require information about the radiative cooling properties of 
the suspended dust. These types of calculations are used in general climate models 
for Mars as well as science and engineering calculations for surface system thermal 
design and scientific measurements on the ground. The optical properties of the dust 
in the infrared are extremely important in any equation of radiative equilibrium on 
Mars, since dust is the only major source of opacity besides CO 2 in the infrared, 
and the CO 2 is mostly limited to a small portion of the infrared spectrum. 

Figure 3.9 clearly represents the improvement in fit of our model over those 
of previous studies [Toon et al. , 1977, Pollack et al. , 1995, Clancy et al. , 1995, 
Forget et al. , 1997]. The quantity plotted is the square of the error between the 
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calculated and actual brightness temperatures, using the optical properties, particle 
size distributions, temperatures, and opacities as specified by each study. The par- 
ticle size distributions assumed in each case are shown in Figure 3.7 and the values 
for some of the other free parameters used in the studies are shown in Table 3.3. 

Toon et al. [1977] concluded that montmorillonite 219b, while it was the best 
fit of the materials tested in the study, could not be the only material in the Mars 
dust, because it was too absorbing in the 500 cm -1 region. The fits to our selected 
spectra using montmorillonite 219b and using the particle size distribution assumed 
by Toon et al. ( r e // = 2.75, u e jj = .4182) are shown in Figure 3.10. Note the severe 
discrepancies caused by the 500 cm -1 absorption of montmorillonite. 

Pollack et al. [1995] sought to obtain a particle size distribution which could 
fit observations in the visible while at the same time still fit the infrared. Using 
the same montmorillonite 219 optical properties, but a new particle smaller size 
distribution (r e // = 2.75, v e jj = .4182) produces the fits shown in Figure 3.11. 

Clancy et al. [1995] proposed still smaller particles and proposed the material 
palagonite [Roush, 1994] as the analog for the optical properties. He found that the 
fit in the 500 cm -1 region was much better than Toon and Pollack, but that the fit 
to the major 1075 cm -1 feature was degraded. Brightness temperature using these 
assumptions are given in Figure 3.12. 

Finally, Forget et al. [1997] empirically modified the optical properties of mont- 
morillonite by setting ni=constant=.5 and nr=constant=2.2 in the whole region 
from 200-550 cm -1 . These optical properties produce a better fit in that wavenum- 
ber region, but they are not physically consistent with the Kramers Kronig relations. 
Fits produced using his optical properties are shown in Figure 3.13 

Comparing our fits (see Figure 3.1) with all of the above demonstrates the 
spectral improvement in the fit to the data using the iterated optical properties, 
surface temperatures, and opacities, and using the particle size distribution of Pol- 
lack et al. . 

The second category of calculations for which improved optical constants are 
necessary involve atmospheric heating rate or flux calculations. These quantities 
are hemispherical mean quantities, so it is not immediately obvious what the effect 
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Figure 3.9: This figure clearly illustrates the improvement in fit to the 
Mariner 9 brightness temperature spectra by comparing results calcu- 
lated using indices of refraction and particle size distributions from var- 
ious previous studies. Here we plot the square of the error between 
calculated and actual brightness temperatures at each wavelength. The 
total error integrated over the spectrum is given by Chi and is indicated 
for each study. 
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Figure 3.10: Fits to the selected brightness temperature spectra gener- 
ated using the particle size distribution, optical properties, and temper- 
atures assumed by Toon et al. [1977]. 
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Figure 3.11: Fits to the selected brightness temperature spectra gener- 
ated using the particle size distribution, optical properties, and temper- 
atures assumed by Pollack et al. [1995]. 
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Figure 3.12: Fits to the selected brightness temperature spectra gener- 
ated using the particle size distribution, optical properties, and temper- 
atures assumed by Clancy et al. [1995]. 
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Figure 3.13: Fits to the selected brightness temperature spectra gener- 
ated using the particle size distribution, optical properties, and temper- 
atures assumed by Forget et al. [1997]. 
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of slightly altering optical properties will be on the integrated fluxes. As a demon- 
stration of the effects of our improved optical properties, we have calculated both 
upward and downward fluxes in the atmosphere for each orbit by integrating the 
upward and downward intensities from the multilayer radiative transfer code over 
all vew angles and over all wavelengths. In Figure 3.14, we compare the fluxes re- 
sulting from our optical properties (solid lines) to those of montmorillonite (dotted 
lines). It is immediately obvious that the effect is large, with montmorillonite un- 
derpredicting the downward flux at the surface by nearly half in the dustier orbits 
and overpredicts the upward flux at the top of the atmosphere by nearly ten per- 
cent. Figure 3.15 shows the comparison over all orbits of the upward and downward 
fluxes. 
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Figure 3.14: Upward and downward fluxes computed using current mont- 
morillonite optical properties (solid lines) compared with those computed 
using the new optical properties presented in this thesis for dusty condi- 
tions such as those found during Mariner 9’s orbit 8. The current model 
underpredicts the downwelling radiation at the surface by nearly half, 
and overpredicts the upwelling radiation at the top of the atmosphere by 
nearly ten percent. 
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Figure 3.15: Upward and downward fluxes computed using current mont- 
morillonite optical properties (stars) compared with those computed us- 
ing the new optical properties presented in this thesis (X’s) for all selected 
orbits. The current model underpredictis the downwelling radiation at 
the surface by up to nearly 50radiation at the top of the atmosphere by 
up to nearly ten percent. 
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Table 3.1: Orbit data computed using the Pollack et al. [1995] size dis- 
tribution 


Spec 

r 

V 


T t0 p, K 

^bot? K 

^surfj K 

mm 

a 

1 

1.339 

0.500 

0.417 

214.177 

243.746 

271.877 

■HI 

0.642 

2 

1.292 

0.500 

0.904 

199.792 

259.488 

260.417 

0.660 

0.642 

3 

0.904 

0.500 

0.969 

196.730 

258.875 

265.273 

0.660 

0.642 

4 

0.739 

0.500 

0.387 

191.079 


282.057 

0.660 

0.642 

5 

0.292 

0.500 

0.996 

179.595 

249.974 

276.924 

0.660 

0.642 

6 

0.270 

0.500 

0.982 

178.460 



0.660 

0.642 


0.320 

0.500 

0.527 

174.716 

245.766 

282.262 

0.660 

0.642 


0.270 

0.500 

0.505 

168.652 

241.217 

279.393 

0.660 

0.642 


0.208 

0.500 

0.980 

168.652 

241.217 


0.660 

0.642 

10 

0.085 

0.500 

0.974 

169.684 



0.660 

0.642 


Table 3.2: Sensitivity of model to parameters, given by rms error (Kelvin) 
in brightness temperature over all wavenumbers and orbits. 


Par am 

Ave. 

008 

056 


080 

138 

146 


174H 

174L 

200 

Baseline 

■ 

2.31 

1.15 

1.17 

1.43 

0.76 

0.82 

1.23 

1.08 

0.94 

1.23 

T +10% 


1.47 

1.51 

1.48 

3.27 

0.86 

0.85 

2.12 

1.36 

0.85 

1.13 

T -10% 

1 

3.82 

1.23 

1.24 

1.50 

0.91 

0.99 

0.86 

1.54 

1.17 

1.35 

Ti +5K 

1.98 

5.34 

2.46 

1.85 

1.61 

0.89 

0.99 

0.81 

1.32 

1.08 

1.30 

Ti -5 K 

1.66 

2.23 

2.73 

2.10 

3.23 

0.86 

0.84 

1.87 

1.21 

0.87 

1.17 

T a +5K 

4.33 

4.31 

3.27 

3.98 

2.50 

5.21 

5.34 

4.03 

5.33 

5.83 

6.39 

Tg- 5K 

4.11 

2.04 

3.96 

4.54 

4.59 

5.19 

5.12 

5.66 

4.61 

4.87 

4.67 

a +10% 

1.81 

1.53 

2.07 

1.89 

4.33 

1.14 

1.12 

2.55 

1.55 

1.03 

1.19 

<7 -10% 

1.95 

5.23 

1.10 

1.14 

2.38 

0.82 

0.90 

0.91 

1.93 

1.15 

1.34 

r m +10% 

1.37 

1.57 

1.66 

1.53 

2.71 

0.92 

0.95 

1.78 

0.96 

0.95 

1.19 

r m -10% 

1.42 

3.74 

0.96 

1.04 

1.18 

0.73 

0.81 

0.81 

1.45 

1.02 

1.28 

nr +10% 

1.49 

2.77 

1.36 

1.22 

2.72 

0.76 

0.86 

1.28 

1.01 

0.93 

1.23 

nr -10% 

1.44 

3.14 

1.38 

1.37 

1.73 

0.88 

0.87 

1.60 

1.39 

1.00 

1.22 

ni +10% 

1.17 

2.03 

1.23 

1.25 

1.84 

0.78 

0.83 

1.45 

1.12 

0.91 

1.20 

ni -10% 

1.18 

2.64 

1.11 

1.13 

1.11 


0.83 

1.02 

1.11 

0.98 

1.26 


















Table 3.3: Parameters used in various studies in the literature 


Study 

Year 

Optical properties 

Teff 

fim 

Veff r m 
fim 

a 

fim 


Toon 

1977 

Montmorillonite 219b 

2.75 

.4182 

1.3 

.548 

Drossart 

1991 

? 

1.24 

.25 

.709 

.472 

Moroz 

1993 

? 

1.6 

.25 

? 

? 

Santee 

1993 

No scattering 

.709 

0.16 

.5 

.387 

Pollack 

1995 

Montmorillonite 219b 

1.85 

0.51 

0.66 

0.642 

Clancy 

1995 

Palagonite 

1.8 

0.8 

.33 

.387 

Forget 

1997 

Mod. montmorillonite 

2.75 

.4182 

1.3 

.548 

Ockert-Bell 

1997 

(visible) 

1.85 

0.51 

0.66 

0.42 

Snook 

1998 

inverted 

1.85 

0.51 

0.66 

0.42 



4. Conclusions 


Our method of directly inverting thermal emission spectra of Mars to determine 
the optical properties, opacity, surface temperatures, and particle size distribution 
of the suspended Mars dust has been demonstrated to be successful in fitting the 
Mariner 9 data to an accuracy that exceeds the accuracy of the data itself. We have 
found that a particle size distribution with an effective radius of reff = 1.85/mi and 
an effective variance of v e jj = 0.51 enables convergence to realistic values for surface 
temperatures (260-281 K) and optical depth (1.34 - 0.18). 
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